home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / LIBRARY / BPL70N16 / ARISOURC.ZIP / FPSIN.ASM < prev    next >
Assembly Source File  |  1993-03-07  |  8KB  |  164 lines
  1.  
  2. ; *******************************************************
  3. ; *                                                     *
  4. ; *     Turbo Pascal Runtime Library Version 7.0        *
  5. ; *     Real Sine/Cosine                                *
  6. ; *                                                     *
  7. ; *     Copyright (C) 1989-1993 Norbert Juffa           *
  8. ; *                                                     *
  9. ; *******************************************************
  10.  
  11.              TITLE   FPSIN
  12.  
  13.  
  14. CODE         SEGMENT BYTE PUBLIC
  15.  
  16.              ASSUME  CS:CODE
  17.  
  18. ; Externals
  19.  
  20.              EXTRN   RealAdd:NEAR,RealSub:NEAR,RealPoly:NEAR,RealMulNoChk:NEAR,
  21.              EXTRN   HaltError:NEAR,RealTrunc:NEAR,RealMul:NEAR,RealFloat:NEAR
  22.              EXTRN   RealReduceMP:NEAR,RFrac:FAR,RealMulFNoChk:NEAR
  23.  
  24. ; Publics
  25.  
  26.              PUBLIC  RSin,RCos
  27.  
  28. ;-------------------------------------------------------------------------------
  29. ; RSin and RCos are different entries into the same routine that computes the
  30. ; sine as well as the cosine, as cosine is the same as the sine with an offset
  31. ; of pi/2. The argument is reduced to the range -pi/4..+pi/4. Depending on the
  32. ; octant into which the argument falls and the function to be computed, the
  33. ; cosine or sine of the reduced argument is computed by polynomial approxi-
  34. ; mations to the sine/cosine. To prevent excessive loss of precision in the
  35. ; argument reduction, the mapping of the argument to the reduced argument is
  36. ; done via a pseudo extended precision calculation. However this extended
  37. ; precision calculation is only available for arguments inside the range
  38. ; -1.68e+9..1.68e+9. Outside this range, a simple argument reduction is used.
  39. ; This causes a sudden drop in accuracy if the absolute value of the argument
  40. ; gets larger than 1.68e+9. For arguments bigger than 6.9e12, there is a total
  41. ; loss of precision. Outside this range, Sin always returns 0, Cos returns 1.
  42. ; The approximating polynomial for the sine on -pi/4..pi/4 is:
  43. ;
  44. ; P(x) := ((((2.476053850842e-8 *x² + 2.755533965768e-6)*x² -
  45. ;             1.984126372865e-4)*x² + 8.333333325083e-3)*x² -
  46. ;             1.666666666663e-1)*x²*x + x
  47. ;
  48. ; This approximation has a theoretical maximum relative error of 1.13e-14.
  49. ; Maximum observed error when evaluated in REAL arithmetic is 1.04e-13.
  50. ;
  51. ; The approximating polynomial for the cosine on -pi/4..pi/4 is:
  52. ;
  53. ; P(x) := ((((-2.715314622037e-7 *x² + 2.479862035332e-5)*x² -
  54. ;              1.388887879411e-3)*x² + 4.166666651281e-2)*x² -
  55. ;              4.999999999923e-1)*x² + 1.000000000000e+0
  56. ;
  57. ; This approximation has a theoretical maximum relative error of 1.41e-13.
  58. ; Maximum observed error when evaluated in REAL arithmetic is 1.49e-12.
  59. ;
  60. ; INPUT:     DX:BX:AX  argument
  61. ;
  62. ; OUTPUT:    DX:BX:AX  sin of argument
  63. ;
  64. ; DESTROYS:  AX,BX,CX,DX,SI,DI,Flags
  65. ;-------------------------------------------------------------------------------
  66.  
  67. RCos         PROC    FAR
  68.              PUSH    BP                ; save caller's frame pointer
  69.              AND     DH, 7Fh           ; abs(x), since cos (x) = cos (abs(x))
  70.              MOV     BP, 2             ; load function code for cosine
  71.              JMP     $cos_entry        ; goto sine/cosine computation
  72. RCos         ENDP
  73.  
  74.              ALIGN   4
  75.  
  76. RSin         PROC    FAR
  77.              PUSH    BP                ; save caller's frame pointer
  78.              XOR     BP, BP            ; load function code for sine
  79. $cos_entry:  MOV     CH, 80h           ; load mask to extract sign bit
  80.              AND     CH, DH            ; extract sign bit
  81.              XOR     DH, CH            ; x := abs(x)
  82.              PUSH    CX                ; save sign
  83. $value_ok:   MOV     CX, 04E81h        ; load
  84.              MOV     SI, 0836Eh        ;  real constant
  85.              MOV     DI, 022F9h        ;   4/pi
  86.              CMP     AL, 9Fh           ; x >= 2^30 ?
  87.              JAE     $big_val          ; x/(pi/4) too big to be stored in long
  88.              PUSH    AX                ; save x
  89.              PUSH    BX                ;  on
  90.              PUSH    DX                ;   stack
  91.              CALL    RealMulFNoChk     ; compute x/(pi/4)
  92.              XOR     CH, CH            ; signal truncation desired
  93.              CALL    RealTrunc         ; n = Trunc (x/(pi/4))
  94.              ROR     AL, 1             ; if quadrant odd, set CF=1
  95.              ROL     AL, 1             ; restore register, CF=1 if odd quadrant
  96.              ADC     AX, 0             ; increment quadrant
  97.              ADC     DX, 0             ;  if quadrant odd
  98.              ADD     BP, AX            ; q := quadrant + flag
  99.              CALL    RealFloat         ; convert quadrant to REAL value
  100.              POP     DI                ; get
  101.              POP     SI                ;  back
  102.              POP     CX                ;   x
  103.              TEST    BP, 3             ; determine octants that use cos,sin
  104.              PUSHF                     ; save sin/cos flag
  105.              PUSH    BP                ; save q
  106.              MOV     BP, OFFSET c1     ; pointer to reduction constant table
  107.              CALL    RealReduceMP      ; compute x + n*(pi/2) to full precision
  108.              MOV     CX, 5             ; use five coefficients
  109.              XOR     SI, SI            ; signal P(x²)*x+x wanted
  110.              POP     BP                ; get back q
  111.              POPF                      ; get back sin/cos flag
  112.              JPE     $sine             ; octants 0, 3, 4, 7 take sine
  113. $cosine:     MOV     DI,OFFSET CS_COEFF; pointer to 1st coeff. of cos approx.
  114.              INC     SI                ; signal P(x²) wanted
  115.              CALL    RealPoly          ; P(x²) in DX:BX:AX
  116.              MOV     CX, 00081h        ; load
  117.              XOR     SI, SI            ;  floating-point
  118.              MOV     DI, SI            ;   constant 1.0
  119.              CALL    RealAdd           ; 1 + P(x²)
  120.              JMP     $make_sign        ; put appropriate sign on result
  121. $sine:       MOV     DI,OFFSET SN_COEFF; pointer to 1st coeff. of sin approx.
  122.              CALL    RealPoly          ; P(x²)*x+x in DX:BX:AX
  123. $make_sign:  POP     CX                ; get sign mask
  124.              INC     BP                ; determine if
  125.              TEST    BP, 4             ;  result has to negated
  126.              JZ      $ende             ; result ok
  127.              XOR     CH, 80h           ; negate sign flag
  128. $ende:       POP     BP                ; restore caller's frame pointer
  129.              OR      AL, AL            ; result zero ?
  130.              JZ      $end_sin          ; yes, done
  131.              XOR     DH, CH            ; make sign bit
  132. $end_sin:    RET                       ; done
  133.  
  134. $big_val:    SUB     CL, 3             ; 1/(2*pi) in DI:SI:CX
  135.              CALL    RealMulNoChk      ; x/(2*pi)
  136.              CALL    RFrac             ; frac (x/(2*pi))
  137.              MOV     CX, 02183h        ; load
  138.              MOV     SI, 0DAA2h        ;  constant
  139.              MOV     DI, 0490Fh        ;   2*pi
  140.              CALL    RealMulNoChk      ; 2*pi * frac (x/(2*pi)) = remainder
  141.              JMP     $value_ok         ; reduction to 0..2*pi done
  142.  
  143. SN_COEFF     DB      067h,               0B1h,0D4h  ; -2.476053850842e-8
  144.              DB      06Eh,05Ch,08Bh,0B6h,0EBh,038h  ;  2.755533965768e-6
  145.              DB      074h,0B5h,09Ch,0FCh,00Ch,0D0h  ; -1.984126372865e-4
  146.              DB      07Ah,044h,086h,088h,088h,008h  ;  8.333333325083e-3
  147.              DB      07Eh,0A9h,0AAh,0AAh,0AAh,0AAh  ; -1.666666666663e-1
  148. CS_COEFF     DB      06Bh,               0C7h,091h  ; -2.715314622037e-7
  149.              DB      071h,034h,0B7h,0A1h,006h,050h  ;  2.479862035332e-5
  150.              DB      077h,02Eh,00Ah,058h,00Bh,0B6h  ; -1.388887879411e-3
  151.              DB      07Ch,018h,0A0h,0AAh,0AAh,02Ah  ;  4.166666651281e-2
  152.              DB      07Fh,0EFh,0FFh,0FFh,0FFh,0FFh  ; -4.999999999923e-1
  153.  
  154. C1           DB      080h,000h,000h,000h,00Fh,0C9h  ; pi/4-eps =
  155. C2           DB      070h,0c2h,068h,021h,0a2h,0dah  ;      eps =
  156.  
  157. RSIN         ENDP
  158.  
  159.              ALIGN   4
  160.  
  161. CODE         ENDS
  162.  
  163.              END
  164.